#include <algorithm/analytical_model/CommonPara.h>
#include <algorithm/analytical_model/Constants.h>
#include <algorithm/analytical_model/HeliostatPara.h>
#include <algorithm/analytical_model/HeliostatWithGlassPara.h>
#include <algorithm/analytical_model/ScenePara.h>
#include <algorithm/analytical_model/SimulationPara.h>
#include <algorithm/analytical_model/algorithm/ModelPara.h>
#include <algorithm/analytical_model/algorithm/PeakFluxParam.h>
#include <algorithm/analytical_model/algorithm/geometry/AllGeometry.h>
#include <algorithm/analytical_model/algorithm/geometry/Utils.h>
#include <algorithm/analytical_model/convolution/HuangPara.h>
#include <algorithm/shared/data/scene/Heliostat.h>
#include <algorithm/shared/data/scene/Receiver.h>
#include <algorithm/shared/data/scene/Sun.h>
#include <algorithm/shared/data/topology/geometry/Geometry.h>
#include <algorithm/shared/data/topology/geometry/Rectangle3D.h>
#include <algorithm/shared/data/world/sun/sun_position/SunPosition.h>
#include <utils/HasValue.h>
#include <utils/ParaUtils.h>
#include <utils/config.h>

#include <cmath>
#include <numbers>

using namespace solar;

namespace
{
    void computeHuangSigma(SimulationPara* simulation_para, ScenePara* scene_para,
                           HeliostatPara* heliostat_para, Sun* sun, Heliostat* heliostat,
                           Receiver* receiver, HuangPara* ptr)
    {
        const auto dni = SOLAR_GET_VALUE(ScenePara, DNI_, *scene_para);
        const double sigma_sun = (3.7648 - 3.8413 * (dni / 1000 - 1) +
                                  1.5923 * 100 * (dni / 1000 - 1) * (dni / 1000 - 1)) /
                                 std::numbers::sqrt2_v<double>;
        const auto sigma_tracking =
            SOLAR_GET_VALUE(HeliostatPara, sigma_tracking_, *heliostat_para);
        const auto sigma_slope_error =
            SOLAR_GET_VALUE(HeliostatPara, sigma_slope_error_, *heliostat_para);
        const auto sun_position_ptr = SOLAR_GET_VALUE(Sun, sun_position_, *sun);
        const auto azimuth = SOLAR_GET_VALUE(SunPosition, azimuth_, *sun_position_ptr);
        const auto elevation = SOLAR_GET_VALUE(SunPosition, elevation_, *sun_position_ptr);
        const auto sun_direction_inv = angleToDirection(azimuth, elevation);
        const auto geometry_heliostat = SOLAR_GET_VALUE(Heliostat, geometry_, *heliostat);
        const auto geometry_receiver = SOLAR_GET_VALUE(Receiver, geometry_, *receiver);
        const auto heliostat_normal = ::solar::getNormal(geometry_heliostat.get());
        const auto heliostat_center = ::solar::getCenter(geometry_heliostat.get());
        const auto cos_theta = std::abs(
            solar::dot(sun_direction_inv.get_unit_vector(), heliostat_normal.get_unit_vector()));
        const double sigma_sx_pow_2 = 4 * sigma_slope_error.x() * sigma_slope_error.x() +
                                      (1 - cos_theta) * (1 - cos_theta) / (8 * cos_theta) *
                                          sigma_slope_error.y() * sigma_slope_error.y();
        const double sigma_sy_pow_2 = 4 * sigma_slope_error.y() * sigma_slope_error.y() +
                                      (1 - cos_theta) * (1 - cos_theta) / (8 * cos_theta) *
                                          sigma_slope_error.x() * sigma_slope_error.x();
        const double sigma_tx_pow_2 =
            sigma_sun * sigma_sun + sigma_sx_pow_2 + 4 * sigma_tracking.x() * sigma_tracking.x();
        const double sigma_ty_pow_2 =
            sigma_sun * sigma_sun + sigma_sy_pow_2 + 4 * sigma_tracking.y() * sigma_tracking.y();
        const auto focus_center =
            ::solar::getFocusCenter(geometry_heliostat.get(), geometry_receiver.get());
        const auto reflectivity = SOLAR_GET_VALUE(HeliostatPara, reflectivity_, *heliostat_para);
        const auto distance_heliostat_focus_center_pow_2 =
            (heliostat_center - focus_center).squared_length();
        const double peak_flux_param =
            peakFluxParam(ModelType::kConvolution, focus_center, heliostat_center, heliostat_normal,
                          0, dni, reflectivity) /
            distance_heliostat_focus_center_pow_2 / ((1 - cos_theta) * (1 - cos_theta));

        SOLAR_SET_VALUE(HuangPara, sigma_x_,
                        std::sqrt(sigma_tx_pow_2 / distance_heliostat_focus_center_pow_2), *ptr)
        SOLAR_SET_VALUE(HuangPara, sigma_y_,
                        std::sqrt(sigma_ty_pow_2 / distance_heliostat_focus_center_pow_2), *ptr)
        SOLAR_SET_VALUE(HuangPara, flux_param_, peak_flux_param, *ptr)
    }

    void computeHuangSigmaWithGlass(SimulationPara* simulation_para, ScenePara* scene_para,
                                    HeliostatPara* heliostat_para, Sun* sun, Heliostat* heliostat,
                                    Receiver* receiver, HuangPara* ptr)
    {
        const auto dni = SOLAR_GET_VALUE(ScenePara, DNI_, *scene_para);
        const double sigma_sun = (3.7648 - 3.8413 * (dni / 1000 - 1) +
                                  1.5923 * 100 * (dni / 1000 - 1) * (dni / 1000 - 1)) /
                                 std::numbers::sqrt2_v<double>;
        const auto sigma_tracking =
            SOLAR_GET_VALUE(HeliostatPara, sigma_tracking_, *heliostat_para);
        const auto sigma_glass_down_slope_error =
            SOLAR_GET_VALUE(HeliostatPara, sigma_slope_error_, *heliostat_para);
        const auto sun_position_ptr = SOLAR_GET_VALUE(Sun, sun_position_, *sun);
        const auto azimuth = SOLAR_GET_VALUE(SunPosition, azimuth_, *sun_position_ptr.get());
        const auto elevation = SOLAR_GET_VALUE(SunPosition, elevation_, *sun_position_ptr.get());
        const auto sun_direction_inv = angleToDirection(azimuth, elevation);
        const auto geometry_heliostat = SOLAR_GET_VALUE(Heliostat, geometry_, *heliostat);
        const auto geometry_receiver = SOLAR_GET_VALUE(Receiver, geometry_, *receiver);
        const auto heliostat_normal = ::solar::getNormal(geometry_heliostat.get());
        const auto heliostat_center = ::solar::getCenter(geometry_heliostat.get());
        const auto cos_theta = std::abs(
            solar::dot(sun_direction_inv.get_unit_vector(), heliostat_normal.get_unit_vector()));
        const auto cos_theta_pow_2 = cos_theta * cos_theta;
        const auto sin_theta_pow_2 = 1 - cos_theta_pow_2;
        const auto heliostat_glass_para_ptr =
            SOLAR_GET_VALUE(HeliostatPara, heliostat_glass_para_, *heliostat_para);
        const auto refractivity =
            SOLAR_GET_VALUE(HeliostatWithGlassPara, refractivity_, *heliostat_glass_para_ptr);
        const auto refractivity_pow_2 = refractivity * refractivity;
        const auto sigma_glass_up_slope_error = SOLAR_GET_VALUE(
            HeliostatWithGlassPara, sigma_glass_up_slope_error_, *heliostat_glass_para_ptr);
        const auto sqrt_common = std::sqrt(std::abs(refractivity_pow_2 - sin_theta_pow_2));
        const auto sigma_glass_up_slope_error_x_pow2 =
            sigma_glass_up_slope_error.x() * sigma_glass_up_slope_error.x();
        const auto sigma_glass_up_slope_error_y_pow2 =
            sigma_glass_up_slope_error.y() * sigma_glass_up_slope_error.y();
        const auto sigma_glass_down_slope_error_x_pow2 =
            sigma_glass_down_slope_error.x() * sigma_glass_down_slope_error.x();
        const auto sigma_glass_down_slope_error_y_pow2 =
            sigma_glass_down_slope_error.y() * sigma_glass_down_slope_error.y();

        // Compute sigma_sx_avg

        const double sigma_sx_pow_2_part1_sqrt = 1 - cos_theta / sqrt_common;
        const double sigma_sx_pow_2_part1 = sigma_sx_pow_2_part1_sqrt * sigma_sx_pow_2_part1_sqrt;
        const double sigma_sx_pow_2_part2 =
            sigma_glass_up_slope_error_x_pow2 +
            ((2 * refractivity_pow_2 - sin_theta_pow_2) / (2 * refractivity * sqrt_common) - 1) *
                sigma_glass_up_slope_error_y_pow2;
        const double sigma_sx_pow_2_part3 = 4 * sigma_glass_down_slope_error_x_pow2;
        const double sigma_sx_pow_2_part4 =
            ((4 * refractivity_pow_2 - 2 * sin_theta_pow_2) / (refractivity * sqrt_common) - 1) *
            sigma_glass_down_slope_error_y_pow2;
        const double sigma_sx_pow_2_part5_sqrt = (1 - sqrt_common) / cos_theta;
        const double sigma_sx_pow_2_part5 = sigma_sx_pow_2_part5_sqrt * sigma_sx_pow_2_part5_sqrt;
        const double sigma_sx_pow_2_part6 =
            sigma_glass_up_slope_error_x_pow2 +
            (1 - cos_theta) * (1 - cos_theta) / (2 * cos_theta) * sigma_glass_up_slope_error_y_pow2;
        const double sigma_sx_pow_2 = sigma_sx_pow_2_part1 * sigma_sx_pow_2_part2 +
                                      sigma_sx_pow_2_part3 + sigma_sx_pow_2_part4 +
                                      sigma_sx_pow_2_part5 * sigma_sx_pow_2_part6;

        // Compute sigma_sy_avg

        const double sigma_sy_pow_2_part1_sqrt = 1 - cos_theta / sqrt_common;
        const double sigma_sy_pow_2_part1 = sigma_sy_pow_2_part1_sqrt * sigma_sy_pow_2_part1_sqrt;
        const double sigma_sy_pow_2_part2 =
            sigma_glass_up_slope_error_y_pow2 +
            ((2 * refractivity_pow_2 - sin_theta_pow_2) / (2 * refractivity * sqrt_common) - 1) *
                sigma_glass_up_slope_error_x_pow2;
        const double sigma_sy_pow_2_part3 = 4 * sigma_glass_down_slope_error_y_pow2;
        const double sigma_sy_pow_2_part4 =
            ((4 * refractivity_pow_2 - 2 * sin_theta_pow_2) / (refractivity * sqrt_common) - 1) *
            sigma_glass_down_slope_error_x_pow2;
        const double sigma_sy_pow_2_part5_sqrt = (1 - sqrt_common) / cos_theta;
        const double sigma_sy_pow_2_part5 = sigma_sy_pow_2_part5_sqrt * sigma_sy_pow_2_part5_sqrt;
        const double sigma_sy_pow_2_part6 =
            sigma_glass_up_slope_error_y_pow2 +
            (1 - cos_theta) * (1 - cos_theta) / (2 * cos_theta) * sigma_glass_up_slope_error_x_pow2;
        const double sigma_sy_pow_2 = sigma_sy_pow_2_part1 * sigma_sy_pow_2_part2 +
                                      sigma_sy_pow_2_part3 + sigma_sy_pow_2_part4 +
                                      sigma_sy_pow_2_part5 * sigma_sy_pow_2_part6;

        const double sigma_tx_pow_2 =
            sigma_sun * sigma_sun + sigma_sx_pow_2 + 4 * sigma_tracking.x() * sigma_tracking.x();
        const double sigma_ty_pow_2 =
            sigma_sun * sigma_sun + sigma_sy_pow_2 + 4 * sigma_tracking.y() * sigma_tracking.y();

        const auto focus_center =
            ::solar::getFocusCenter(geometry_heliostat.get(), geometry_receiver.get());
        const auto reflectivity = SOLAR_GET_VALUE(HeliostatPara, reflectivity_, *heliostat_para);
        const auto distance_heliostat_focus_center_pow_2 =
            (heliostat_center - focus_center).squared_length();
        const double peak_flux_param =
            peakFluxParam(ModelType::kConvolution, focus_center, heliostat_center, heliostat_normal,
                          0, dni, reflectivity) /
            distance_heliostat_focus_center_pow_2 / ((1 - cos_theta) * (1 - cos_theta));

        SOLAR_SET_VALUE(HuangPara, sigma_x_,
                        std::sqrt(sigma_tx_pow_2 / distance_heliostat_focus_center_pow_2), *ptr)
        SOLAR_SET_VALUE(HuangPara, sigma_y_,
                        std::sqrt(sigma_ty_pow_2 / distance_heliostat_focus_center_pow_2), *ptr)
        SOLAR_SET_VALUE(HuangPara, flux_param_, peak_flux_param, *ptr)
    }

} // namespace

template <>
void solar::computeSigma(SimulationPara* simulation_para, ScenePara* scene_para,
                         HeliostatPara* heliostat_para, CommonPara* common_para, Sun* sun,
                         Heliostat* heliostat, Receiver* receiver, HuangPara* ptr)
{
    const auto heliostat_type = SOLAR_GET_VALUE(HeliostatPara, type_, *heliostat_para);
    if (heliostat_type == HeliostatType::kWithGlass)
    {
        computeHuangSigmaWithGlass(simulation_para, scene_para, heliostat_para, sun, heliostat,
                                   receiver, ptr);
    }
    if (heliostat_type == HeliostatType::kWithoutGlass)
    {
        computeHuangSigma(simulation_para, scene_para, heliostat_para, sun, heliostat, receiver,
                          ptr);
    }
}
